library(tidyverse) #ggplot2, dplyr, tibble, etc.
library(plotly) #use ggplotly(*plotName*) for interactive plots with scoll-over IDs
library(knitr)
library(ggfortify)
library(Rtsne)
library(usethis) #git integration
library(glmmTMB)
library(lme4)
library(ggpubr) #ggarrange() figure wraps
library(plyr)
library(ggeffects) #ggpredict
library(rstanarm) #stan models
library(shinystan) #stan model evaluation >launch_shinystan_demo()
library(loo) #loo() to compare fits between bayesian models
library(MuMIn) #dredge, model averaging
library(RColorBrewer) ; brewer.pal(n=8, name = "Dark2")#; display.brewer.pal(n = 8, name = "Dark2")
## [1] "#1B9E77" "#D95F02" "#7570B3" "#E7298A" "#66A61E" "#E6AB02" "#A6761D"
## [8] "#666666"
theme_set(theme_classic())

Create full dataframe (df).

df <- data.frame(read.csv("Nestlings1989-2018_no2011_WITH_AdultMeasurements_2001-2015.csv", h=T))
summary(df[,c(2:23)])
##       year       population          species           lifeStage        
##  Min.   :1989   Length:4372        Length:4372        Length:4372       
##  1st Qu.:1995   Class :character   Class :character   Class :character  
##  Median :2002   Mode  :character   Mode  :character   Mode  :character  
##  Mean   :2002                                                           
##  3rd Qu.:2009                                                           
##  Max.   :2018                                                           
##                                                                         
##       id                 age               nest             nestName        
##  Length:4372        Min.   :    0.00   Length:4372        Length:4372       
##  Class :character   1st Qu.:   10.00   Class :character   Class :character  
##  Mode  :character   Median :   24.00   Mode  :character   Mode  :character  
##                     Mean   :   28.77                                        
##                     3rd Qu.:   27.00                                        
##                     Max.   :22014.00                                        
##                     NA's   :336                                             
##   sexSummary            tarsus        sevenPrimary       tailLength    
##  Length:4372        Min.   :  0.00   Min.   :    0.0   Min.   :  0.00  
##  Class :character   1st Qu.: 43.80   1st Qu.:   27.0   1st Qu.:  7.00  
##  Mode  :character   Median : 57.00   Median :   97.0   Median : 47.00  
##                     Mean   : 49.43   Mean   :   92.5   Mean   : 49.82  
##                     3rd Qu.: 59.90   3rd Qu.:  121.0   3rd Qu.: 69.00  
##                     Max.   :144.00   Max.   :14550.0   Max.   :195.00  
##                     NA's   :30       NA's   :100       NA's   :123     
##     calcAge             billNT             billD           billW       
##  Min.   :    0.00   Min.   : -0.1224   Min.   : 0.00   Min.   :  0.00  
##  1st Qu.:   13.00   1st Qu.: 16.1000   1st Qu.:11.10   1st Qu.: 11.60  
##  Median :   24.98   Median : 21.4000   Median :12.40   Median : 12.60  
##  Mean   :   29.91   Mean   : 19.8106   Mean   :11.95   Mean   : 12.14  
##  3rd Qu.:   26.31   3rd Qu.: 23.4000   3rd Qu.:13.20   3rd Qu.: 13.40  
##  Max.   :29344.26   Max.   :158.0000   Max.   :20.40   Max.   :112.60  
##  NA's   :551        NA's   :48         NA's   :60      NA's   :60      
##       TEC             head            skull          wingChord    
##  Min.   : 0.00   Min.   :  0.00   Min.   :-52.80   Min.   :  0.0  
##  1st Qu.:24.90   1st Qu.: 60.70   1st Qu.: 35.60   1st Qu.:103.0  
##  Median :33.20   Median : 73.90   Median : 40.30   Median :157.0  
##  Mean   :30.29   Mean   : 67.64   Mean   : 37.33   Mean   :142.2  
##  3rd Qu.:36.00   3rd Qu.: 78.10   3rd Qu.: 42.50   3rd Qu.:181.0  
##  Max.   :58.50   Max.   :101.70   Max.   : 93.90   Max.   :329.0  
##  NA's   :53      NA's   :59       NA's   :55       NA's   :356    
##   exp7Primary          weight     
##  Min.   :   0.00   Min.   :  0.0  
##  1st Qu.:  39.00   1st Qu.:200.0  
##  Median :  56.00   Median :345.0  
##  Mean   :  58.26   Mean   :291.7  
##  3rd Qu.:  73.00   3rd Qu.:400.0  
##  Max.   :2038.00   Max.   :615.0  
##  NA's   :2673      NA's   :120

Add a new column with the number of occurrences of each individual in the dataset.

df$indivCounts <- as.numeric(ave(df$id, df$id, FUN = length))
#sorted by ascending year and lifeStage (add sorting conditions with: ', variableName')
df <- df %>% arrange(year, lifeStage)
#make lifeStage a factor
df$lifeStage <- as.factor(df$lifeStage)
#check n for both life stages
length(df$lifeStage[df$lifeStage == "nestling"]) #4086 nestlings
## [1] 4086
length(df$lifeStage[df$lifeStage == "adult"]) #286 adults
## [1] 286
#create column with numbered observations for reference
df$obs <- seq(1,length(df$nestlingMeasurementID))
#make sex a factor
df$sexSummary <- as.factor(df$sexSummary)
levels(df$sexSummary)
## [1] "F" "M"

Create dataframe subsets.

#only individuals with multiple measurements (mm) 
mm <- subset(df, df$indivCounts>1)
#only nestlings (ne)
ne <- subset(df, df$lifeStage == "nestling")
#only adults (ad)
ad <- subset(df, df$lifeStage == "adult")
#subset morphometric data from .ne
nestlMorph.ne <- ne[,c(1:23)] #n=4086
#drop NAs
nestlMorph.ne <- drop_na(nestlMorph.ne)#n=1413
nestlMorph.ne$sexSummary <- as.factor(nestlMorph.ne$sexSummary)
class(nestlMorph.ne$sexSummary)
## [1] "factor"
levels(nestlMorph.ne$sexSummary)
## [1] "F" "M"
#filter ages (also reduced by NA removal for unknown sexes)
nestlMorph.ne <- nestlMorph.ne %>% filter(between(age, 23, 27))
#create nestling bill morphometric dataframe
nestlBills.ne <- nestlMorph.ne[,c(15:18)]
tarsusByWeight <- nestlMorph.ne %>%  
  ggplot(aes(x=weight,y=tarsus,color=sexSummary,label=id))+
  geom_point()
ggplotly(tarsusByWeight)
skullByWeight <- nestlMorph.ne %>%  
  ggplot(aes(x=weight, y=skull, color = sexSummary, label = id))+
  geom_point()
ggplotly(skullByWeight)
billPCA <- prcomp(nestlBills.ne)
summary(billPCA)
## Importance of components:
##                           PC1    PC2     PC3    PC4
## Standard deviation     2.5344 0.9077 0.76513 0.5385
## Proportion of Variance 0.7908 0.1014 0.07207 0.0357
## Cumulative Proportion  0.7908 0.8922 0.96430 1.0000
billPCA
## Standard deviations (1, .., p=4):
## [1] 2.5344278 0.9077234 0.7651260 0.5384864
## 
## Rotation (n x k) = (4 x 4):
##               PC1         PC2         PC3        PC4
## billNT -0.5308950  0.08039039 -0.82350303 -0.1831138
## billD  -0.1795232  0.30231519 -0.06245138  0.9340647
## billW  -0.1414367  0.90873294  0.24765403 -0.3047419
## TEC    -0.8160377 -0.27631026  0.50656694 -0.0335406
billPCA.ne <- tibble(PC1=billPCA$x[,1],
                     PC2=billPCA$x[,2],
                     sex=nestlMorph.ne$sexSummary,
                     id=nestlMorph.ne$id)

billPCA_plot <- billPCA.ne %>%
  ggplot(aes(x=PC1,y=PC2,color=sex,label=id))+
  geom_point()
ggplotly(billPCA_plot)
morphPCA <- prcomp(nestlMorph.ne[,11:23])

morphPCA.plot <- tibble(PC1=morphPCA$x[,1],
                        PC2=morphPCA$x[,2],
                        sex=nestlMorph.ne$sexSummary,
                        id=nestlMorph.ne$id)
summary(morphPCA)
## Importance of components:
##                            PC1     PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     53.7309 26.2520 5.45909 4.19985 2.75474 2.28349 1.85734
## Proportion of Variance  0.7922  0.1891 0.00818 0.00484 0.00208 0.00143 0.00095
## Cumulative Proportion   0.7922  0.9813 0.98949 0.99433 0.99642 0.99785 0.99879
##                            PC8     PC9    PC10    PC11    PC12     PC13
## Standard deviation     1.64247 0.87822 0.66903 0.50423 0.47736 2.94e-15
## Proportion of Variance 0.00074 0.00021 0.00012 0.00007 0.00006 0.00e+00
## Cumulative Proportion  0.99953 0.99974 0.99987 0.99994 1.00000 1.00e+00
morphPCA
## Standard deviations (1, .., p=13):
##  [1] 5.373088e+01 2.625199e+01 5.459091e+00 4.199851e+00 2.754744e+00
##  [6] 2.283486e+00 1.857336e+00 1.642469e+00 8.782189e-01 6.690293e-01
## [11] 5.042252e-01 4.773575e-01 2.940383e-15
## 
## Rotation (n x k) = (13 x 13):
##                       PC1           PC2          PC3          PC4          PC5
## tarsus       -0.038179895 -0.0041571604  0.005423045  0.120275476  0.366072576
## sevenPrimary -0.180750280  0.5007468830  0.142762776  0.311802527 -0.718103471
## tailLength   -0.169665443  0.3659809253  0.555711993 -0.711826105  0.121406957
## calcAge      -0.009755084  0.0292277133  0.010712123 -0.013512652 -0.014560321
## billNT       -0.018667171  0.0191472815  0.005238650  0.020614391  0.086587176
## billD        -0.009043039  0.0009808141 -0.005119253  0.008156555 -0.007620977
## billW        -0.005325301 -0.0001540695  0.006722487  0.015456294  0.002797935
## TEC          -0.024821963  0.0234086209 -0.003332407  0.012766991  0.107772601
## head         -0.045274211  0.0406639533 -0.027252951  0.096631723  0.205755069
## skull        -0.020452248  0.0172553324 -0.023920545  0.083864732  0.097982468
## wingChord    -0.212187158  0.4938662426  0.197633543  0.511596188  0.506516463
## exp7Primary  -0.131362064  0.4958120660 -0.791512479 -0.319731594  0.077758697
## weight       -0.933363781 -0.3490539194 -0.060687124 -0.014515091 -0.040588008
##                        PC6          PC7          PC8           PC9
## tarsus       -0.3839684220  0.139908741  0.823167452 -0.0593058494
## sevenPrimary -0.2437940532  0.020134803  0.150395150  0.0096874574
## tailLength   -0.0561074066  0.054020382  0.004834509 -0.0073901265
## calcAge       0.0006989063 -0.018584139  0.007881656  0.0025541398
## billNT       -0.2503866342 -0.179780817 -0.166879870  0.0002973483
## billD        -0.0288488520 -0.066656717 -0.048722349 -0.2476934611
## billW        -0.0739067474 -0.060008910 -0.086528860 -0.9562300860
## TEC          -0.3451388317 -0.674472290 -0.078289287  0.1003405420
## head         -0.6229716734  0.016993897 -0.382415366  0.1014748445
## skull        -0.2778328417  0.691466187 -0.304126079  0.0011343025
## wingChord     0.3663868389 -0.054623473 -0.128056173 -0.0002542232
## exp7Primary   0.0318983036  0.008833085  0.025970132 -0.0146684138
## weight        0.0365215872 -0.001528300 -0.006702304  0.0042229533
##                       PC10         PC11         PC12          PC13
## tarsus       -0.0373501140  0.011609262  0.019037700  1.238924e-16
## sevenPrimary -0.0098638877  0.014317600 -0.021248184  1.119139e-16
## tailLength    0.0117819828  0.025549436 -0.009160179 -7.667362e-17
## calcAge       0.0132875073 -0.755419223  0.653675739  6.533814e-16
## billNT       -0.9222981215 -0.095524280 -0.093342817 -7.585426e-17
## billD        -0.1091091636  0.626340942  0.725642728  6.659476e-16
## billW         0.0837301651 -0.163407019 -0.187191368  2.443571e-16
## TEC           0.2514661116  0.006843275 -0.014155675 -5.773503e-01
## head          0.2554900619  0.009773461  0.015998322  5.773503e-01
## skull         0.0040239503  0.002930186  0.030153998 -5.773503e-01
## wingChord     0.0089653382  0.003504284 -0.003161622 -4.357466e-17
## exp7Primary  -0.0059954773  0.005145221 -0.010006872  7.870149e-17
## weight       -0.0001802609 -0.005463276 -0.004859557 -2.439204e-17
indMorphPCA <- cbind(morphPCA$x[,1], nestlMorph.ne$id)

morphPCA <- morphPCA.plot %>%
  ggplot(aes(x=PC1,y=PC2,color=sex,label=id))+
  geom_point()
ggplotly(morphPCA)

Comparisons of bill measurements by sex

Including nestlings 20-30 days in age as per calcAge (n=538)

BillNT

nestlMorph.ne %>% 
  ggplot(aes(x=sexSummary,y=billNT))+
  geom_boxplot()

t.test(nestlMorph.ne$billNT[nestlMorph.ne$sexSummary=="M"], nestlMorph.ne$billNT[nestlMorph.ne$sexSummary=="F"])
## 
##  Welch Two Sample t-test
## 
## data:  nestlMorph.ne$billNT[nestlMorph.ne$sexSummary == "M"] and nestlMorph.ne$billNT[nestlMorph.ne$sexSummary == "F"]
## t = 3.7865, df = 318.26, p-value = 0.0001826
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.2874067 0.9091245
## sample estimates:
## mean of x mean of y 
##  22.09316  21.49490

BillW

nestlMorph.ne %>% 
  ggplot(aes(x=sexSummary,y=billW))+
  geom_boxplot()

t.test(nestlMorph.ne$billW[nestlMorph.ne$sexSummary=="M"], nestlMorph.ne$billW[nestlMorph.ne$sexSummary=="F"])
## 
##  Welch Two Sample t-test
## 
## data:  nestlMorph.ne$billW[nestlMorph.ne$sexSummary == "M"] and nestlMorph.ne$billW[nestlMorph.ne$sexSummary == "F"]
## t = 5.4804, df = 353.99, p-value = 8.076e-08
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.3285849 0.6964170
## sample estimates:
## mean of x mean of y 
##  12.89538  12.38288

BillD

nestlMorph.ne %>% 
  ggplot(aes(x=sexSummary,y=billD))+
  geom_boxplot()

t.test(nestlMorph.ne$billD[nestlMorph.ne$sexSummary=="M"], nestlMorph.ne$billD[nestlMorph.ne$sexSummary=="F"])
## 
##  Welch Two Sample t-test
## 
## data:  nestlMorph.ne$billD[nestlMorph.ne$sexSummary == "M"] and nestlMorph.ne$billD[nestlMorph.ne$sexSummary == "F"]
## t = 4.3053, df = 332.26, p-value = 2.197e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.1790865 0.4804187
## sample estimates:
## mean of x mean of y 
##  12.70203  12.37227

TEC

nestlMorph.ne %>% 
  ggplot(aes(x=sexSummary,y=TEC))+
  geom_boxplot()

t.test(nestlMorph.ne$TEC[nestlMorph.ne$sexSummary=="M"], nestlMorph.ne$TEC[nestlMorph.ne$sexSummary=="F"])
## 
##  Welch Two Sample t-test
## 
## data:  nestlMorph.ne$TEC[nestlMorph.ne$sexSummary == "M"] and nestlMorph.ne$TEC[nestlMorph.ne$sexSummary == "F"]
## t = 4.9907, df = 313.39, p-value = 9.994e-07
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.6710488 1.5445208
## sample estimates:
## mean of x mean of y 
##  34.53278  33.42500

BillW

df %>% 
  ggplot(aes(x=lifeStage,y=tarsus,color=sexSummary))+
  geom_boxplot()
## Warning: Removed 30 rows containing non-finite values (stat_boxplot).

df %>% 
  ggplot(aes(x=lifeStage,y=billW,color=sexSummary))+
  geom_boxplot()
## Warning: Removed 60 rows containing non-finite values (stat_boxplot).

df %>% 
  ggplot(aes(x=lifeStage,y=billD,color=sexSummary))+
  geom_boxplot()
## Warning: Removed 60 rows containing non-finite values (stat_boxplot).

df %>% 
  ggplot(aes(x=lifeStage,y=billNT,color=sexSummary))+
  geom_boxplot()
## Warning: Removed 48 rows containing non-finite values (stat_boxplot).

#create morphometric matrix
nestlMorph.mat <- as.matrix(nestlMorph.ne[,-c(1:10,22)])
#run tSNE
nestlMorph.tsne <- Rtsne(nestlMorph.mat)

nestlMorphTsne.plot <- data.frame(x=nestlMorph.tsne$Y[,1],
                                 y=nestlMorph.tsne$Y[,2],
                                 nest=nestlMorph.ne$nest,
                                 population=nestlMorph.ne$population,
                                 sex=nestlMorph.ne$sexSummary,
                                 year=nestlMorph.ne$year)

Clustering by year

ggplot(nestlMorphTsne.plot, aes(x=x,y=y,color=as.factor(year)))+
  geom_point(size = 2)+
  xlab("Dimension 1")+
  ylab("Dimension 2")+
  theme(#legend.position="none",
        legend.title = element_text(size = 18, face = "bold"),
        axis.title.x = element_text(size = 18, face = "bold"),
        axis.text = element_text(size = 16),
        axis.title.y = element_text(size = 18, face = "bold"))

Clustering by sex

Clustering by nest ID

Relationships between calculated ages and morphometrics

CalcAge by weight

cAge_morphs.ne <- ne %>% filter(between(calcAge, 1,50))
cAge_morphs.ne$year <- as.factor(cAge_morphs.ne$year)

#filter out weights <=1
cAge_morphs.filter.ne <- cAge_morphs.ne %>% filter(weight >=1)
summary(cAge_morphs.ne$calcAge)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00   13.00   24.98   20.11   26.28   50.00
summary(cAge_morphs.filter.ne$weight)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     5.0   240.0   350.0   304.3   397.0   500.0
#ageWtCheck <- as.data.frame(cbind(Wt = cAge_morphs.ne$weight, cAge = cAge_morphs.ne$calcAge))
cAge_morphs.filter.ne %>% 
  ggplot(aes(x = calcAge, y = weight, color = as.factor(year)))+
  geom_point(shape = 1)+
  theme(legend.position = "none")

Zoom in on cluster in the 24-30 day range

wtMidCluster <- cAge_morphs.filter.ne %>% filter(between(calcAge, 24,30))
wtMidCluster %>% 
  ggplot(aes(x = calcAge, y = weight, color = as.factor(year)))+
  geom_point(shape = 1)+
  theme(legend.position = "none")

Decimal points come from the formula for age calculation. Should we be rounding to nearest day?

CalcAge by tail length

cAge_morphs.ne %>% 
  ggplot(aes(x = calcAge, y = tailLength, color = as.factor(year)))+
  geom_point(shape = 1)+
  theme(legend.position = "none")
## Warning: Removed 90 rows containing missing values (geom_point).

Zoom in on cluster in the 24-30 day range

tailMidCluster <- cAge_morphs.ne %>% filter(between(calcAge, 24,30))
tailMidCluster %>% 
  ggplot(aes(x = calcAge, y = tailLength, color = year))+
  geom_point(shape = 1)+
  theme(legend.position = "none")
## Warning: Removed 8 rows containing missing values (geom_point).

CalcAge by 7th primary

cAge_morphs.ne %>% 
  ggplot(aes(x = calcAge, y = sevenPrimary, color = as.factor(year)))+
  geom_point(shape = 1)+
  theme(legend.position = "none")
## Warning: Removed 74 rows containing missing values (geom_point).

Zoom in on cluster in the 24-30 day range

primMidCluster <- cAge_morphs.ne %>% filter(between(calcAge, 24,30))
primMidCluster %>% 
  ggplot(aes(x = calcAge, y = sevenPrimary, color = year))+
  geom_point(shape = 1)+
  theme(legend.position = "none")
## Warning: Removed 1 rows containing missing values (geom_point).

CalcAge by tarsus

cAge_morphs.ne %>% 
  ggplot(aes(x = calcAge, y = tarsus, color = as.factor(year)))+
  geom_point(shape = 1)+
  theme(legend.position = "none")
## Warning: Removed 13 rows containing missing values (geom_point).

Zoom in on cluster in the 24-30 day range

tarsMidCluster <- cAge_morphs.ne %>% filter(between(calcAge, 24,30))
tarsMidCluster.filter <- tarsMidCluster %>% filter(between(tarsus, 20,100))
tarsMidCluster.filter %>% 
  ggplot(aes(x = calcAge, y = tarsus, color = year))+
  geom_point(shape = 1)+
  theme(legend.position = "none")